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1. PROTECT SUMMARY 


PROTECT OVERVIEW: 

The contract "Near Field Interaction of Microwave Signals with a Bounded Plasma 
Plume" (under Cooperative Agreement NCC3-127) was awarded to the University of 
Texas at Austin in January 1989. The objective of the project is to study the effect of the 
arcjet thruster plume on the performance of an on-board satellite reflector antenna. The 
total amount awarded, $50,000 + $50,000 = $100,000, covers the two-year period from 
January 20, 1989 to January 19, 1991. 


Technical Monitor: Dr. Afroz Zaman 

M.S. 54-8 

NASA Lewis Research Center 
21000 Brookpark Road 
Cleveland, OH 44135 
Tel: (216)-433-3415 

FAX: (2 1 6)-433-637 1 


Principal Investigators: Dr. Hao Ling and Dr. Gary A. Hallock 

Department of Electrical and Computer Engineering 
The University of Texas at Austin 
Austin, TX 78712-1084 
Tel: (5 1 2)-47 1-1710 (H.L.) 

(5 12)-47 1-4965 (G.A.H.) 

FAX: (5 12)-47 1-5445 


Other Investigators: Mr. Hyeongdong Kim and Mr. Bjorn Birkner 

Graduate Students 


This Final report covers the two-year period from Janurary 20, 1989 to January 19, 
1991, with special emphasis on the work performed during the second year. 
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TRAVEL 

(1) Drs. Ling and Hallock visited NASA Lewis Research Center on February 3, 1989 for 
a kick-off meeting on the project. They each presented an one-hour talk on the plasma 
modeling and antenna aspects of the project. 

(2) Dr. Zaman visited the University of Texas at Austin on August 25, 1989 for an interim 
review of the project. She also toured the microwave and plasma fusion facilities. 

(3) Drs. Ling and Hallock visited NASA Lewis Research Center on February 16, 1990 to 
report the progress and discuss future directions of the project. 

(4) Drs. Ling and Hallock visited NASA Lewis Research Center on July 30, 1990 to 
report on the project and discuss the experimental aspects with Mr. Bjorn Birkner, 
who participated in the summer graduate student program. 
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FIRST- YEAR ACCOMPLISHMENTS 


(1) Plasma modeling has been carried out. Electron density data supplied by NASA on 
the arcjet plume has been used for an in-depth modeling of the spatial plume profile as 
a function of arcjet power and mixture ratio. Spatial permittivity distribution of the 
plume based on the cold plasma model has been derived for the microwave 
propagation modeling. 

(2) The finite temperature effect in the arcjet plume have been studied and found to be 
negligible for the range of operating conditions encountered in laboratory arcjets. 

(3) A two-dimensional ray-tracing code for an arbitrary inhomogeneous plasma has been 
developed. It will calculate the ray path as well as the phase, amplitude and 
polarization of the associated ray optics field. 

(4) The numerical results from the 2-D ray tracing algorithm have been compared to exact 
solutions available for several special media, e.g.. Maxwell's fisheye and a 
quadratically varying ionosphere. Good agreement is obtained given enough precision 
in the numerical calculations. 

(5) A two-dimensional reflector pattern computation code has been implemented. A ray- 
tube integration scheme is used to calculate the far field from the ray optics field 
transmitted through the plume. The far field performance degradation of reflectors in 
the presence of the plume has been calculated in 2-D. It is found that in addition to 
gain loss, beam broadening and sidelobe deterioration, a significant beam squint can 
occur due to refraction of ray paths in the plume. 

(6) In order to determine the lower frequency bound of validity of the ray approach, the 
scattering of a plane wave by an inhomogeneous circular cylinder was formulated 
using ray optics. The ray results were compared to (i) a special medium where the 
exact solution to this problem exists, and (ii) a layered cylinder model where again the 
exact solution can be found. The agreement between the ray results and the exact 
solutions is excellent when the cylinder radius is above 3 free-space wavelengths. 


4 


(7) A three-dimensional ray-tracing code has being developed. The 3-D calculation of 
field amplitude is significantly more complicated than the 2-D version. In addition, the 
rotation of the polarization vector due to the non-planar nature of the ray path is 
properly taken into account Preliminary results of the 3-D implementation have been 
verified by the exact results available in Maxwell's fisheye. 

(8) A three-dimensional reflector pattern computation code has been implemented. A ray- 
tube integration scheme is used to calculate the far field from the ray optics field 
transmitted through the plume. The far field performance degradation of reflectors in 
the presence of the plume has been calculated in 3-D. Reflector performance 
degradation such as beam squint, beam broadening, sidelobe deterioration and cross 
polarization can be evaluated quantitatively. 

(9) Plans for possible experimental validation of the numerical modeling has been 
initiated. Both near field and far field measurement setups are under investigation. 
We are also looking into using microwave interferometer for non-invasive plume 
diagnostics. 
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SECOND-YEAR ACCOMPLISHMENTS 


(1) Collisional effects in the arcjet plasma which can result in wave attenuation were 
investigated. They were found to be negligible for the range of operating conditions 
encountered in the laboratory arcjet plume. 

(2) Using the ray modeling codes developed in the first year, 3-D reflector performance 
degradation information such as gain loss, pointing error, sidelobe levels and cross 
polarization have been tabulated for different plume-to-reflector configurations. It was 
found that the plume effect on the reflector performance is strongly dependent on the 
position of the reflector with respect to the arcjet thruster. 

(3) Plans for experimental validation of the numerical ray modeling were initiated at the 
end of the first year. It was decided that reflector far-field measurements would be 
diffcult to realize in the laboratory, especially in view of the vacuum environment 
necessary for arcjet operation. A near- field experiment has therefore been further 
pursued. Since the plume effect is expected to be small, wall reflection inside the 
vacuum facility is an important consideration. Vacuum compatible absorbers were 
considered for damping out wall reflections. (Several sheets of Eccosorb FDS were 
purchased for testing.) In addition, the software time-gating feature of the HP8510 
network analyzer was extensively investigated. By making measurements over a 
broad frequency range and Fourier transforming to the time domain, software gating 
can be applied in the time domain to eliminate the reflection signals. Measurements of 
a dielectric (paraffin) lens were made in the laboratory and the results demonstrated 
that this concept can be applied with good accuracy. 

(4) Our graduate student, Bjorn Birkner, participated in the NASA summer program. He 
interacted closely with both the propulsion group and the antenna group and utilized 
the in-house facility at NASA to set up microwave experiments similar to those 
performed at Texas. The feasibility of the plume/reflector near-field experiment using 
NASA's vacuum facility was further investigated. Measurement accuracy and aliasing 
errors from the software time-gating concept were also considered in more detail. At 
the end of the summer, he assisted Dr. Zaman in putting together a detailed cost 
estimate of carrying out the experiment. 
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(5) Further development in the ray modeling code has been carried out, with the aim of 
supporting the experimental effort. First, the 3-D near-field calculation has been 
implemented, in order to simulated the experimental setup. Next, the variable step size 
feature in the ray tracing routines was implemented to improve the efficiency of the ray 
code. This capability would allow larger reflectors to be modeled. Finally, the effect 
of a sharp plume boundary has been investigated in order to (i) simulate experiments 
carried out in a non-vacuum environment, such as measuring the scattering pattern 
from a dielectric lens, and (ii) model sharp plume boundaries which may occur in high 
power arcjets. Preliminary results showed good agreement when compared with exact 
solutions for simple scatterers. 

(6) In parallel with the experimental effort, an independent validation of the 3-D ray 
tracing code has been carried out. The scattering from an inhomogeneous sphere was 
considered. The series expansion solution to the scattering from a multi-layered 
sphere was implemented to approximate the exact solution to the problem. The 
agreement between the ray results and the reference solution is quite good when the 
sphere diameter is above 6 wavelengths, especially in the forward scattering region. 
This confirms the validity of the ray calculation and further supports our earlier results 
on the reflector performance prediction. 
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CONCLUSION S AND RECOMMENDATIONS 


(1) From measurement data and theoretical calculations made on the arcjet plasma, we find that 
higher order effects such as anisotropy due to the earth's magnetic field, finite plasma 
temperature and attenuation due to collisions are small for low-power arcjets (~ 1 kW). A 
simple collisionless cold plasma model should be adequate in describing the plume plasma. 

(2) Electromagnetic modeling using the approximate ray theory has been verified against various 
exact solutions for both two- and three-dimensional scatterers. We find that the ray theory is 
accurate as long as the size of the plume scatterer is on the order of 5 wavelengths or greater. 
The accuracy of the ray theory is particularly good in the forward scattering region (i.e., the 
transmitted field through the plume). In the reflector antenna study, we are mainly interested 
in this region. 

(3) The far-field patterns of the reflector in the presence of the plume have been computed for 
different reflector-to-plume configurations and the resulting performance degradation have 
been tabulated. As the arcjet power is increased, a gradual degradation of the reflector pattern 
is observed. The most significant and interesting effect is a beam squinting phenomenon 
where the antenna mainbeam tends to tilt away from the plume region. Associated with this 
phenomenon are antenna gain loss and pointing error. The amount of degradation in the 
antenna performance is also found to be strongly dependent on the specific plume-to-reflector 
configuration. Finally, for low-power arcjets in the 1-kW class, the plume effect on the 
antenna performance is almost negligible. 

(4) Experimental verification of the results predicted by the ray theory may be difficult and costly. 
A near-field measurement performed inside a vacuum facility is probably the most feasible 
option, as oppose to carrying out a far field test. Since an anechoic vacuum chamber is 
difficult to achieve, additional signal processing techniques such as time-gating are needed to 
prevent wall reflections from contaminating the measurement. 

(5) The computer modeling capabilities developed in this work can be applied to predict the plume 
impact on microwave links in more realistic antenna and arcjet configurations, such as those 
found on the AT&T TELSTAR-4 communication satellite. 
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2 Plasma Modeling 


As a medium for the propagation of electromagnetic waves, a plasma is refrac- 
tive, lossy, dispersive, nonlinear, and inhomogeneous (and additionally resonant, 
anisotropic, and non-reciprocal if magnetized). 

2.1 Wave Propagat ion 

Wave propagation is governed by Maxwell’s equations. Combining Ampere s and 
Faraday’s laws yields the wave equation 

d — dE 

v x (v X E) + -QtiVoJ + e oVo-fo) = 0 


Assuming: 

1. Linearity (wave fields sufficiently small so current is a linear function of electric 
field) 

2. Medium is homogeneous in space and time 
can Fourier analyze the fields and currents 

£(*,<) = J E{k,u})S i£ - ut) d 3 kduj 

and write an Ohm’s law expression for each Fourier mode 

J(k,u)) = cr(A:,u>) • E(k, u>) 
where a is the conductivity tensor. 

The plasma does not satisfy the spatial uniformity assumption. However, if 
gradients are sufficiently small then locally the above analysis can be applied. An 
eikonel or geometrical optics approximation can be used. 

For a single Fourier mode the wave equation yields 

k x (k x E) + ■ E — e 0 /xo iuE) = 0 

(kk - k 2 I- fwV) E=0 

where I is the unit dyad and l the dielectric tensor 



The plasma properties may be equivalently specified by e or a. For a nontrivial 
solution for E 

det(kk — k 2 I + = 0 

which is the dispersion relation, relating k and w for different waves. Its solution 
also gives the wave polarization. 

If the plasma is not immersed in a magnetic field it can usually be treated as 
an isotropic medium. Then 


a = al 


6 = el 


where <r and e are the scalar conductivity and permittivity. The determinant reduces 
to 

det{kk - k 2 I + u>Voc) = {-k 2 + u>Voe) Woe) 
which is equal to zero if 


o> 2 /i 0 e — k 2 = 0 


or 


u; 2 /x 0 e = 0 

The first equation represents a solution with transverse polarization of the electric 
field (k ■ E = 0) while the second is a longitudinal solution (k x E = 0). 

The transverse solution 

U) 

k = -- s/Zr 
c 

can be expressed in terms of the index of refraction of the medium (u p = u >/k) 

c _ 

N = — = 

v p 


and is of interest here. 

The longitudinal wave dispersion relation is e = 0, which can have a nontrivial 
solution in a plasma. 
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2.2 Plasm a Analy sis 

We consider an infinite, uniform, unmagnetized Lorentz plasma. The Lorentz model 
assumes 

1. Electrons interact with each other only through collective spacecharge forces. 

2. Ions and neutrals can be regarded as a continuous stationary fluid through 
which the electrons move with viscous friction 

We analyze the motion of individual representative particles, ignoring the statistical 

correlations in the electrons position and velocity. 

In addition, we ignore thermal particle motion (T e — T, = 0), which is valid 
when v th < c and v p ~ c (cold plasma). We therefore assume the only particle 
motion comes from the wave fields. This is governed by the Lorentz equation. 

2.2.1 Collisionless. Model 
The electron motion is given by 

m^j- = q{E + vxB) 
at 

The wave fields are related by 

E __ HIo 
H * V e 

which gives 

Fe _ qE _ Vp c 

Fb qvB v v 

Therefore the wave’s magnetic force on the plasma can be neglected, and the Lorentz 
equation reduces to 

dv - 
m— = qF 
dt 

In the Fourier analysis v is harmonic 

—iuimv — qE 



Lorn 


The current density is given by 
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J — qnv 


J = 


e hn~ 

u>m e 


which gives the plasma conductivity ( J — crE) 


e 2 m 

<x = 

uim e 


and permittivity 

i 

e = eo H & 

UJ 


e = e 0 - 


e 2 n 

a > 2 m c 


Defining the plasma frequency 


2 2 2 2 
U>_ = U>. + (jJ: ~ U) = 

P Co"l e 


yields the cold plasma permittivity 




- e o(! - .. 2 ) 


u; 46 


and the dispersion relation 


u; / u* 


k = -Ml 




2.2.2 Collisional Model 

The influence of discrete positive ions and neutral molecules in a plasma can be rep- 
resented to good approximation by including a viscous damping term, proportional 
to velocity, in the electron equation of motion 

dv - 

m = —eE - vmv 
dt 

With Fourier analysis this yields 
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— € 


V — 


m{ v — iu>) 


E 


Using J = -env = crE gives 


e 2 n 


<j — 


m(u — iu>) 


The permittivity is given by 

i 

6 = 60 + — O’ 
U> 


which gives 


e r = 1 


a;; 


u ){iv + a?) 


or 


Cp = (1 - 


iVZ 


V 2 + U)' 


r) + i( 




lt)(v 2 + u> 2 ) 


) 


Now using 


u; 


k = — s/T t = f3 + ia 


gives the collisional dispersion relation 


k = 


u> 


1 - 


u>- 


+ i 


VU)t 


c ^ v 2 + u 2 U)(v 2 + a? 2 ) 


Three frequency regimes can be identified: 

1. Low frequencies u; < u. The high plasma conductivity for DC and low fre- 
quency fields is obtained. 

2. Intermediate frequencies u < u) < u p . Here the plasma will not propagate an 
EM wave, similar to a waveguide in cutoff. The wave penetration is the order 
of the wavelength of a free space wave at the plasma frequency. 

3. High frequencies u; > u> p . The plasma is a relatively low-loss dielectric. This 
is the regime of interest. 
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For the high frequency regime, using 




U) 


\J\ + X ~ 



yields 


and therefore 


k ~ 



ivu: 


2u; 3 (l - a;* /a; 2 ), 


f3 ~ 





or 


a ~ 


a? 

c 



Mp 

2uj 3 


The real refractive index is quite insensitive to collisional damping. The attenuation 
is small except near the plasma frequency. 


2.2.3 Collisions 

For a partially ionized gas 


total ^en t" ^ ci 

where v en is the electron-neutral collision frequency and v ei is the electron-ion col- 
lision frequency. 

The degree of ionization above which Coulomb collisions dominate is strong y 
temperature dependent. At T e ~ 0.1 eV the critical ionization is small, the order of 
10~ 5 — 10 -6 . The arcjet operates with n e ~ 5 x 10 _4 n„ and T e from 0.1 to 0.3 e\ 

(Carney NASA TM 102284). 

For a Maxwellian velocity distribution 

— 8.4 x 10 [T c ^ 3/2 s 

with n in cm- 3 and T e in eV. For n e = 10 10 cm" 3 and T e = 0.1 eV 

i/ ei = 2.7 x 10 7 s - 1 
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The collision frequency for the scattering of charged particles of species a by 
neutrals is 

[ kT, * _i 

V an — U-n Cr c m ^ ^ 

V TUa 

where n n is the neutral density and a an is the cross section, typically ~ 5 x 10 16 cm 2 
and weakly dependent on temperature. For electrons 

i/ en = 2.1 x 10 - 7 n nS fT e a" 1 

with n n in cm" 3 and T e in eV. For n n = 2 x 10 13 cm" 3 and T e = 0.3 eV 

u en = 2.3 x 10 8 a -1 

These calculations are consistent with Carney’s mean free path estimates 

A et = 3 — 200 cm 


A en = 20 - 4000 cm 


Using 


v — ncrv 


v 

A 


and 


ISkTe 

7T m e 


v(T' = 0.1 eV) = 2.1 x 10 7 cm/a 


v(T e = 0.3 eV) = 3.7 x 10 7 cm/a 


gives 

1.1 x 10 5 < v ei < 7 x 10 6 a -1 
9.3 x 10 3 < v m < 1.9 x 10 6 a -1 
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2.2.4 Attenuation 


Using the following parameters 

n = 10 10 cm -3 =*> u) p = 5.65 x 10 e s _1 


1 / = 3 x 10 7 s" 1 



yields, for ui — 1.2tci p , attenuation and skin depth 

a = 0.0628 m _1 


For = 1.5u> p the results are 



15.9 m 


a = 0.0298 m -1 


<5 = — = 33.6 m 

a 

and for u> = 5u^ p 

a = 0.00203 m" 1 

6 = — = 492 m 
a 

Conclusions: 

1. The attenuation is small for the low power arcjet, even close to the plasma 
frequency. 

2. Attenuation could be a factor in signal degradation in higher power arcjets if 
v increases about an order of magnitude. This will depend on how the plasma 
density and temperature scales. 

3. These results can be incorporated into our ray-tracing codes. A perturbation 
technique is most appropriate, where the ray trajectory is unchanged (the 
collisional 0 is unchanged from the collisionless model) but the attenuation a 
is carried along giving exponential attenuation of the ray amplitude. 
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2.3 Plume Elect ron Density Model 


1. We have modeled the plume profile as a function of space and arcjet power, 
based on Carney’s experimental measurements and motivated by NASA mod- 
eling. 

2. We have initially used a simple model allowing linear optimization: n e = 
Aiexp(—a0)/r 2 . 

3. We have also evaluated Reynold’s model using nonlinear optimization: n e = 
A 0 exp{ — [A(l — cos(0)] ;v }/2 2 . 

4. Future work will include modeling other functions to the experimental data, 
to eliminate the on-axis discontinuity and test the ray-tracing code sensitivity 
to the model details. 

5. Our modeling programs allow optimization of almost any function with the 
experimental data available. The code provides the least squared error for 
comparison. 

6. We have also modeled the arcjet power dependence, based on the experimental 
data: n e = n e {P = 1025 W)K P exp{/3 P P). 

7. The work has been done on a DEC microVAX workstation, using IMSL math 
routines and IDL processing/plotting software. 

8. The output of this modeling is the input to the ray-tracing codes. We have also 
used the microVAX capabilities for interactive 3-d display of the ray-tracing 
code output. 

Figures 1 - 5 illustrate this work. The analysis is described in detail in Publi- 
cation 8 (B. W. Birkner M.S. Thesis). 
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2.4 Prop osed Experimental Program 

1. Experimental verification of our present results is a logical next step in our 
program. This should be done in parallel with configuration studies and ad- 
vanced modeling/code development. 

2. Experimental verification is very important because the geometrical optics 
assumptions which underlie the ray-tracing methods may not be well satisfied 
by the plasma plume. 

3. In addition, it appears that the equipment and experimental setup used for 
code validation could be used as a microwave interferometer for plume electron 
density measurements. 

2.4.1 Npnr-FielH Measurements 

1. The issue for satellite communications is the plume effect on the far-field 
antenna radiation pattern. 

2. Far-field measurements are difficult to simulate in the laboratory, however, 

with the plume vacuum requirements. A far-field experiment mig t e possi e 

in the 19 m NASA, tank, but would require dedicated use of the tank an 
would only allow limited angular measurements. 

3. Therefore a near-field experiment seems more appropriate. It will require 
much less space, and allow measurements over the entire antenna pattern. 

4. A near field experiment is also a very sensitive test of our code. If our near- 
field calculations are verified correct then our far-field calculations are very 

accurate. 

5. We propose a near field experiment be undertaken during the summer of 1990. 
One of our graduate students, Bjorn Birkner, will work on-site at NASA Lewis. 
We would hope equipment costs can be kept low by making use of microwave 
components available at NASA Lewis and UT Austin. 

2.4.2 Electron Density Measurements 

1 Experimental measurements of the plasma electron density are important for 
optimizing arcjet performance. The profile is also a fundamental parameter 
in our communication degradation studies. 

2 It appears that the experimental apparatus required for code verification can 
also be used for microwave interferometry. Additional components required 

are: 
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(a) Two microwave lenses to focus the antenna radiation and provide spatial 
resolution. 

(b) Microwave hardware (waveguide, beamsplitters, etc.) to configure the 
system as a Mach-Zehnder interferometer. 

(c) A mounting structure (or beam steering capabilities) to allow positioning 
the beam with respect to the plasma. 

2.4.3 Comparison of Interferometry with Probe Meas urements 

ADVANTAGES: 

1. It is completely non-invasive. 

2. It can be applied to higher power arcjets (and in fact will work better). Lang- 
muir probes will be restricted to the edge region of high power thrusters. 

3. It can provide measurements close to the thruster nozzle (Langmuir probes 
appear to be limited to about 12 cm). Data near the arcjet nozzle are needed 
to understand plume anomalies (n e increases with increased flow rate; expect 
more recombination, etc.). 

4. It might be capable of resolving probe measurement uncertainties, due to 
unknown factors such as the “wetted” current collection area. 

5. Cross-section or even full 3-D profiles can be made in a short time. 

6. It has a fast time response (typically sub-millisecond) for transient measure- 
ments. 

7. It has high noise immunity (excellent environmental isolation). The system 
should operate well in pulsed ignition experiments, for example. 

DISADVANTAGES: 

1. The measurement is inherently chord-integrated and requires Abel or tomo- 
graphic inversion. Langmuir probes provide a highly localized measurement. 

2. Application of interferometry to low power (~ 1 K XV) arcjets appears feasible, 
but nontrivial. It is near the lower limit of plasma size x density needed for 
the measurement. However, all aspects should improve as the arcjet program 
evolves to higher power thrusters. 

3. Interferometry is more complex and costly than probe measurements. 
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2.4.4 Proposed Arcjet Interferometer 
CONSTRAINTS: 

1. To avoid dissipative and nonlinear effects the operating frequency should be 
several times the plasma frequency. 

2. The probing wavelength should be small compared to the plasma diameter, 
to minimize refractive effects. 

3. The probing beam should be small compared to the plasma size, to provide 
good spatial resolution. 

4. The operating frequency must not be too far above the plasma frequency or 
the phase shift will be too small to measure. 

TNTTTAL DESIGN: 

1. 10 GHz operating frequency (The plasma frequency is ~ 1 GHz for n e = 
10 10 cm- 3 . 

2. Obtain a phase shift of ~ 10 degrees for a 20 cm chord-integrated average 
density of 10 10 cm -3 . 

3. Spatial resolution will be the order of 3 cm, the beam wavelength. This is 15 
% of a 20 cm plume diameter. 

4. A heterodyne configuration is being evaluated. 

A block diagram of the proposed interferometer is shown in Fig. 6. 
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Figure 6 Microwave Interferometer for Plasma Density Measurements. 
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3. ELECTROMAGNETIC MODELING 


3.1. New Developments in Ray Code 

Emphasis : Support experimental effort 

1 . 3-D Near Field Calculation 

2. Modeling Scatterer Boundary 

3. Variable Step Size Ray Tracing 
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1. 3-D Near Field Calculation 


f = 10 GHz 
A1 = 1E13 



100 cm 


Approach: 

• Launch a uniformly spaced set of rays from source antenna. 

• Obtain a set of non-uniformly spaced exit rays on the X ' out. 

• Use "near-field" ray tube integration to calculate the field distribution 

onlouL 

Near-field ray tube integration: 

E(r) =| dS' { i(0|Jg(r,r') [nxH(r')] + VVg(f,?>[nxH(r')] 

I k 2 

J Ray tube 

+Vg(?,r')x[nxE(?')] ) 
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3-D Results 


Data Sampling: 


r 4J-F-70 


y- 



r=U 


Reflector antenna 
projection onto Xout 


Zout aperture 


Data sampling at 4> = 0 degree : 
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2. Modeling Scatterer Boundary 



Sharp discontinuity 
in refractive index 


Motivation: Comparison with preliminary experiment 

Approach : 

1) Propagation through continuous inhomogeneous media 



2) Reflection/ refraction at boundary discontinuity 
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3. Variable Step Size Ray Tracing 


Fixed step size in ray tracing algorithm 

Too much computation time !!! 



Variable step size 

=^> Use large step size when refractive index 
variation is small 




Algorithm 



r h (s 0 +h) 

?h/2( s 0 + h) 


Local truncation error - D n = | fh(So+h) ■ iWCso+h) 

r h (s 0 +h) : ray trajectory using step size h 

Th/ 2 (so+h) : ray trajectory using step size h/2 
ERRLIM : error tolerance 


Based on local error, 

if Dn < ERRLIM, r h (s 0 +h) is accepted and increase stepsize 
if Dn > ERRLIM, reduce step size and try again 

Details in Conte and Boor, Elementary Numerical Analysis, 3rd ed., 
p 366-372,1981 
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3.2. Validation of 3-D Ray Results 


1) Layered model solution 

Problem statement: Find the scattered field of a spherically inhomogeneous 
sphere by modelling it as a multi-layered sphere. 


Details in J. R. Wait, "Electromagnetic scattering from a radially 
inhomogeneous sphere," Appli. Sci. Res., section B, vol. 10, pp. 441-450, 1963 



Ek = EoeJ kz 


Far-field scattered field: 


E* =^eJ kr cos(t»P(0) 

0 kr 

E! = e ^ kr sirKp Q(B) 
0 kr 


where 


P(0)=X 

n=l 


(2n+l) 

n(n+l) 


[B n -^P n (cos0) + C n 
d0 


P„(cos0) 
sin 0 


] 


Q<9) = X 

n=l 


(ii+l) B P B (cosO ) + qjLp^cos 0) ] 
n ( n+1 ) sin 0 d0 
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J p (kai) - jApJ D (kai) 

Hn X (kai)-jAoHi, 2), (kai) 

J B (fcai)-jSbJn(fcai) 

H? X (kai)-j5bH ( n 2/ (kai) 


k = coVeoPo 
Ap = ZjAlo 


5b = Vfilo 


Input impedance and admittance from the outside can be calculated from the 
following recurrence relationship* 


y n _ k 

"m-1 


[Jntkam-ijNpQa m) - NnQam-ijJnOam)] +^Z^[N n (to m -i0 n (kfl m ) - Jn(ka m -i)Nn(katn)] 
jQ)E [J n (ka m . 1 )Np(ka I p) - Np(ka m .i)N„(ka in )] +^Zg,lN„(ka ni -i^p(ka m ) - Jn(Ka m -i)N„(ka I p)] 


tJn(ka I p.i)Np(ka m ) - ^(ka^n^n,)] + ^T“Ym[N n (kani-i)Jp(k a m) - J n (ka m ,)Np(ka in )] 
= k . — “ “ 

^ [J n (ka m -i)N n (kani) - ^(ka^NnOtam)] +^Y&[Np(ka,p.idn(kanO - Jn(ka n ,-i)N n (kan 1 )] 


with m = M, M-l, . 1 


rn j_Jn(kaM) 


/s' 

yn 1 Jn(kaM) 
M jnJn(kaM) 
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2) Ray formulation 


The ray scattering formulation for validation is based on Huygens' principle as in 
two-dimensional case. A dense grid of geometrical optics (GO) rays representing the 
incident plane wave is shot toward the inhomogeneous object. The shot rays are 
traced numerically, considering the discontinuous boundary and propagation inside 
the object. Whenever the rays cross the Huygens' surface S, the contributions of the 
rays to the scattered field are calculated by a ray-tube integration scheme. 



At discontinuous boundary: 

The ray tangential vector, amplitude, polarization and wavefront curvature of 
refflected ( or refracted ) are calculated from those of the incident ray using 
Snell's law, Fresnel's formulations for reflection (or refraction ) coefficient, 
and phase matching. 

Inside the object: 

The ray trajectory, the associated phase, amplitude and polarization of the ray 
field are traced numerically. 
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Formulation of ray tube integration 

: Integrate over ray-tube projection on Huygens' surface 


3-D radiation integral in the far field: 


E s = ^[0Ae + <frA^ 


1 1 

= Jkf e&rl 

471 1 r *^ ube 


- * 

•4> 

xE(r) + r| 

" A * 

0 

xH(r) 

•/projection ' 

on S 

i 

0 

1 

. ♦. 

I 


r = r a + r 

r a = XgX+yay+ZaZ 
r = xx+yy+zz 
v = v x x+v y y+v£ 

n = n x x-m y y-Hi z ?=-^x4^-y4^z 

k = k(cos}> sin0 x+sin<}> sin0 y-tcosB z) 

4» = cos0 cos}) x+cos0 sin<t> y-sin0 z 


0 = -sin<}>x+cos)) y 


E(r) = E(r») exp{ -jkt •(?-?*) } 


H(r) = H(r a ) exp{ -jkt •(?-?,) } = 


VXE^,) 


exp{ -jkt'fr-rn) } 
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'Ae 


= *[ 
4k 1 

J V 


ray tube 
projection 
on S 


eJkTJ 


([ ~ 
4> 


e 


xE(t«)+tJ 


6 

L 5 J 


xHfaje^^aVn ds 


= IM-fJI i 


4k 


6 


xE(r a ) + T| 

/s * 

0 

/*N. 

1 

xH(ra) 



/ 


n (AA) S(0,<t>) 


= i^eS4 B "’ 

4 * [Bj 


(Area)proj S(0,<t>) 


S(0, <t>) = 


1 

(Area)proj 


J raytub* 
projection 
on S 


eJWc-vJ^r'ds' 


Be = ( ^xe + n9xH)«n = [-nxEjCoaJ) - n y E#in$ + n&ysinty + Erfo^)] 

+ TjtrixCHiCoSSsiiKl) + HjSinG) - ^(Hjcose + H x sin0) + n^ycosflcos}) - HxCO^sirnJ))] 


B 0 = (0xE + r$xH)*n 

= n x (E z cos 0sin 0 + E y sin 0) - n y (E z cos 0 + ExSin 0) + n z (E y cos 0cos 4> - E x cos 0sin <|>) 


+ TitnxHzCOS)) + n y H^in4> - n^ysiiKj) + HxCOS<t>)] 


/A (AreaJin, l N 

(Area)proj — ~ zrzi^ ^ 

|DFl 2 n*v 

(Area> n * initial ray-tube size 

S(0, 4 )) = 2 ^ — , u=kA sin [cos’ M(£-v)*n}] 
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3) Comparisons 


Permittivity profile: 


e r (r) = d + b(l-c cos^) 


where 

w 1+E r (0) 

b = "T~ 
_ l-e r (0) 

1+Er(0) 
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Scattering Cross Section (dB X 2 ) S Scattering Cross Section (dB X 2 ) 


(a) 


40 




-90 -60 -30 0 30 60 90 

PHI (DEG) 


Comparion of the scattering cross section of an inhomogeneous sphere of 
diameter 6X calculated using ray optics and that obtained from the multi-layered 
solution, (a) £r(0)=0.5. (b) C((0)«3.0. 
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ARCJET PLASMA PLUME EFFECT ON A 


MICROWAVE REFLECTOR ANTENNA.* 

B.W. Birkner, G.A. Hallock, H. Kim, H. Ling 

The University of Texas at Austin. 

Eighth Topical Conference 
On High Temperature Plasma Diagnostics 
May 6-10 1990 
Poster Presentation F27 


Abstract 

A configuration is described to measure the effect of the plasma pro- 
duced by an arcjet thruster on a reflector antenna. The diagnostic 
consists of a small movable RF probe coupled to a Hewlett-Packard 
HP 8510 network analyzer. Because the perturbation is expected to 
be small (~ 10° phase shift) wall reflections are an important con- 
sideration. To minimize their effect we use Fourier transform gating 
techniques. 

'Supported by the NASA Lewis Research Center 
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INTRODUCTION 


• Electric Propulsion for space applications has regained interest in 
recent years due to availability of large power supplies for satel- 
lites. 

• Arcjet thrusters, one of the simplest devices, can attain specific 
impulse levels above 400 sec. Arcjet thruster propulsion consists 
in heating a flowing gas in a constricted volume by means of an arc 
discharge thus obtaining thrust by heat expansion. See Figure 1. 

• Use of arcjet thrusters for stationkeeping and attitude control can 
potentially extend the lifetime of modern satellites up to ten years. 

• Focus of this work is on the ~ 1 KW class arcjet thruster. 

• Higher power arcjets are being developed for boosting satellites 
from Earth orbit to geosynchronous orbit. 

• The arcjet produces a weakly ionized plasma plume during opera- 
tion. For this technology to be usable on communication satellites, 
the signal degradation to and from the spacecraft due to refrac- 
tion and attenuation of the microwave signals must be within 
acceptable limits. 
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• Establishing arcjet compatibility with microwave communications 
requires both modeling and laboratory experiments, with neither 
approach alone being able to predict the operation of an actual 
system with adequate confidence: 

1. Modeling is limited by plasma parameters which are not ac- 
curately known, and difficulties inherent to wave propagation 
calculations. 

2. Experimental work is limited by laboratory simulation of the 
space environment (boundary conditions, vacuum pressure, 
and wall reflections). 


insulator constricted volume expansion 



Figure 1: Arcjet Thruster 
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• The arcjet produces a weakly ionized plasma which extends 10’s 
of centimeters from the arcjet nozzle. See Figure 2 for arcjet 
geometry and coordinates. 

• The cold plasma model provides a very good approximation to the 
electromagnetic behavior, with index of refraction N given by: 

N = 

where u p is the electron plasma frequency. 

• The plasma frequency dependes on the electron number density ( 
ujp = n e e 2 /m e c 0 ), which has been measured experimentally with 
Langmuir probes. The electron density is ~ 10 10 cm“ 3 close to the 
thruster nozzle, and is well modeled by n e oc exp(—a6)/r 2 for all 
arcjet operating conditions. (See Figure 2). 

• A least squares fit to the data is shown in the contour plot in 
Figure 3 along with the measurement locations. 

• Microwave communication frequencies are typically only a factor 
of 2-20 above the plasma frequency close to the thruster, and the 
plasma behaves as a very dispersive medium. 

• To study the effect of the plasma on the performance of on-board 
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reflector antennas, the arcjet plume is modeled as an inhomoge- 
neous scatterer with a smoothly varying refractive. index. 

• Since the plume size is of the order of several wavelengths, geo- 
metrical optics is used to calculate the wave propagation. 

• Geometrical optics is implemented by means of numerical ray 
tracing. Given the value of the electric and magnetic fields at 
some point, the fields can then be calculated along a ray extend- 
ing from that point. 

• Once ray optics fields are obtained on some near-field aperture of 
the antenna the far-field pattern can be obtained by the Fourier 
Transform. 

• Two-dimensional calculations have been performed indicating the 
plume can cause serious pointing errors, especially for very narrow 
beam antennas. In addition deterioration in both the main beam 
and sidelobes has been observed in the calculations. 
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EXPERIMENT A L VALIDATION 


• Experimental validation is important because the geometrical op- 
tics assumptions which underlie the ray-tracing methods are not 
well satisfied by the plasma plume. For our work plume size ~ A 
whereas geometrical optics requires plume size <C A. 

• Experimental measurement of the effect of the plume on the re- 
flector antenna pattern is planned. 

• Due to size limitations of vacuum chambers a far-field radiation 
pattern is difficult to measure. 

• A near-field radiation measurement is feasible and provides a very 
sensitive test of our code. If near-field calculations are verified 
correct then our far-field calculations are very accurate. 

• The experimental configuration is shown in Figure 4. A 40 cm 
diameter reflector antenna directs collimated radiation at 10 GHz 
through the plasma plume. The diagnostic consisits of a small 
movable RF probe coupled to a Hewlett-Packard HP 8510 network 
analyzer. 

• The arcjet plume will cause the antenna radiation to deflect up- 
ward, since the index of refraction is less than unity within the 
plasma. This changes both the magnitude and phase of the near 
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field distribution, as shown in Figure 5. The effect is small and 
is concentrated towards the region closest to the plume. The in- 
duced phase shift in this region is ~ 8° and the magnitude change 
is of the order of 2.5 decibels. In addition a dark region extending 
for ~2 cm is observed. 

• The arcjet thruster requires a modest vacuum for operation ~ 
10 ~ 4 torr. To reduce the effect of wall reflections we will make use 
of the time domain gating feature of the HP 8510 analyzer. In 
this technique measurements are made over a broad frequency 
range in the AC steady state, and Fourier transformed to the 
time domain via the FFT (Fast Fourier Transform). The FFT is 
applied over a finite bandwidth of positive frequencies. Gating in 
the time domain is then applied to eliminate the reflections, the 
latter being shifted in time by the extra path length traveled. The 
gated data is then transformed back into the frequency domain. 

• Time gating introduces several sources of error: 

1. The finite sampling used for the FFT causes aliasing if the 
extra path length for reflections is too large and/or multiple 
reflections do not die out quickly. 

2. Short path lengths might not be separable in time. 

• Aliasing problem can be reduced by using a sufficiently small 
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chamber so that the first reflections can be gated out, and vac- 
uum compatible microwave absorbing material which can alleviate 
multiple bounce effects. 

• Because the plasma is dispersive the time domain pulses will be 
distorted thus complicating the task of gating out the reflections. 

• A numerical simulation has been done ( shown in Figure 6) where 
the frequency spectrum of the input signal is assumed Gaussian 
with a center frequency of 15 GHz and a FWHM of 2.94 GHz. 

• In Figure 6(b) is shown the magnitude of the transformed pulse 
shape of the input signal. The output of this signal after propa- 
gating through a homogeneous plasma slab of thickness d = 20 cm 
is computed by assuming that each frequency component under- 
goes a phase shift of —{yjd/c)^J 1 — o/p/u; 2 . The results are shown 
in Figure 6(c) for different plasma frequencies. 

• It is shown that as long as the frequency range of measurement is 
sufficiently above the plasma frequency, no significant pulse dis- 
tortion will occur. The observed pulse delay is also not significant. 

• We have set up a test facility consisting of an antenna chamber, 
antennas, and an HP 8510 analyzer. We have verified basic sys- 
tem operation, using time gating on induced reflections. We have 
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successfully compared predicted and measured effects of various 
dielectrics (such as paraffin) to high accuracy. 

• After completing a proof-of-principle work, we plan to set up the 
measurement configuration at NASA Lewis on their arcjet vacuum 
facility. 

• Figure 7 shows the problem of vacuum chamber wall reflections. 
The main beam is illustrated together with two direct reflections - 
the side wall refection and the back wall reflection, and a multiple 
bounce reflection. 

• The following figures illustrate the time gating procedure with 
actual HP 8510 measurements: 

1. Measured frequency response in “clean” environment, Fig- 
ure 8. 

2. Frequency response with a strong reflection, Figure 9. 

3. Frequency response with gated out reflection, Figure 10. 

• Gating errors are of major concern and an illustration of the dis- 
tortion produced on the frequency response by time gating is il- 
lustrated in Figure 11, Figure 12, and Figure 13. Preliminary 
calculations were done to estimate the time gating error in terms 
of the root mean square error over the measuring frequency span: 
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1. For small gating (Figure 12) the frequency response magnitude 
rms-error is ~ 1% , and the phase rms-error is ~ 1.8°. 

2. For large gating (Figure 13) the frequency response magnitude 
rms-error is ~ 0.5% , and the phase rms-error is ~ 0.1°. 

• For a 12-18 GHz reflector feed sweep sampled with 800 points 
we obtain a time domain window of 50 nsec, corresponding to an 
electrical length of 15 m. In addition the time domain 50% pulse 
width is ~ 0.1 nsec, corresponding to an electrical length of ~ 
4 cm. 

• Some considerations on the choice of frequency span are listed 
below: 

1 . Lower frequencies are preferred because more plasma induced 
phase shift can be observed. But for lower frequencies the 
reflector beam is not well as collimated. 

2. Only the interior of the frequency sweep is of concern for the 
plume degradation measurements. 

3. Most of the frequency sweep is used for time gating. 
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▲ x-axis 



Figure 2: Arcjet Thruster Geometry and Coor- 
dinates 
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Figj”™? 3j* Cop fr *ur ( pl ot ( tly» a rc elp'^roif 
density profile. Measurement locations are de- 
noted by x. The contours are given by the least 
squares fit n e = 9.1 x 10 n exp(— O.O190)/r 2 cm" 3 ( 0 in 
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Figure 4: Schematic of the measurement config- 
uration, showing approximate ray trajectories in 
the presence of the plasma. 
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NEAR FIELD MAGNITUDE 




Figure 5: Calculated near-field magnitude and 
phase (referenced to the antenna plane) with and 
without plasma. 
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(a) 


plasma 


zu 

t 

20 cm 


(b) Input Pulse Shape 



(c) Simulated Output Pulse Shape 



Figure 6: Simulated time-domain output of the 
input pulse shape and the output pulse shapes 
after propagation through a "20 cm plasma slab. 




Figure 7: Illustration of Main beam and the dif- 
ferent reflections present in vacuum chamber an- 
tenna measurements. 
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with no reflections. 



60 


with enhanced reflection. 











Frequency Domain 



Figure 11: Simulated FFT with no gating. 
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Arcjet plasma plume effect on a microwave reflector antenna 

B. W. Birkner, G. A. Hallock, H. Kim, and H. Ling 

The University of Texas at Austin, Austin, Texas 78712 
(Presented on 8 May 1990) 

A configuration is described to measure the effect of the plasma produced by an arcjet 
thruster on a reflector antenna. The diagnostic consists of a small movable rf probe coupled 
to a Hewlett-Packard 8510 network analyzer. Because the perturbation is expected to 
be small ( < 10° phVse shift), wall reflections are an important consideration. To minimize 
their effect we use Fourier transform gating techniques. 


I. INTRODUCTION 

The availability of large power supplies on the next 
generation of communication satellites has provided re- 
newed interest in electric propulsion. An arcjet engine is 
one of the simplest devices, and current research indicates 
that specific impulse levels above 400 s are attainable, 
making arcjets a very attractive candidate. The use of 
arcjet thrusters for station keeping and attitude control can 
potentially extend the lifetime of modem communication 
satellites up to ten years. This requires arcjets of the ~ 1 
kW class, which is the focus of the work reported here. 
Higher power arcjets are also being developed for boosting 
satellites from lower Earth orbit to geosynchronous orbit. 
Again the fuel savings over conventional chemical propul- 
sion are very large. 

The arcjet produces a weakly ionized plasma during 
operation. For this technology to be usable on communi- 
cation satellites, the signal degradation to and from the 
spacecraft due to refraction and attenuation of the micro- 
wave signals must be within acceptable limits. Since the 
satellites require totally uninterrupted operation, it is nec- 
essary to understand both steady-state and transient (star- 
tup and shutdown) operation of the thrusters. 

A program of both modeling and laboratory experi- 
ments is required to establish arcjet compatibility with the 
microwave communications. Neither approach alone can 
predict the operation of an actual system with adequate 
confidence. Modeling relies on plasma parameters which 
are not accurately known, and is uncertain because of the 
difficulties inherent in wave propagation calculations. Ex- 
perimental work is limited by the difficulty in laboratory 
simulation of the space environment, particularly bound- 
ary conditions, vacuum pressure, and wall reflections. In 
addition, accurate nonintrusive measurements of the 
plasma parameters are difficult. 

II. PLUME MODELING AND RAY TRACING 

The arcjet produces a weakly ionized plasma plume, 
which extends 10’s of centimeters from the thruster nozzle. 
The cold plasma model provides a very good approxima- 
tion to the electromagnetic behavior, with the index of 
refraction given by N = yf 1 — , where co p is the 
electron plasma frequency. The plasma frequency depends 
on the electron number density (co p = n,£ l /m e e 0 ), which 
has been experimentally measured by Carney using Lang- 


muir probes in a laboratory facility. 2 The electron density 
is — 10 10 cm -3 close to the thruster nozzle, and is well 
modeled by a l/(distance squared) falloff away from the 
nozzle times a factor giving exponentially decreasing den- 
sity with angle away from the thruster axis. A contour 
plot showing the measurement locations and a least- 
squares fit to the data is shown in Fig. 1. The microwave 
communication frequencies are typically only a factor of 
2-20 above the plasma frequency close to the thruster, and 
therefore the plasma behaves as a very dispersive medium. 

To study the effect of the plasma on the performance of 
on-board reflector antennas, the arcjet plume is modeled as 
an inhomogeneous scatterer with a smoothly varying re- 
fractive index. Since the plume size is the order of several 
wavelengths, geometrical optics is used to calculate the 
wave propagation. The trajectory, phase, amplitude, and 
polarization of the high-frequency ray field are computed 
numerically, based on the laws of geometrical optics. 4 

Once the ray optics field transmitted through the 
plume on the near-field aperture is found, the far-field pat- 
tern of the antenna can be calculated by its Fourier trans- 
form. 5 Two-dimensional calculations have been performed, 
and we are now extending our analysis to full 3-D geom- 
etry. Results indicate the plume can cause serious pointing 
errors, especially for very narrow beam antennas. Deteri- 
oration in both the main beam and sidelobes has been ob- 
served in the calculations. 6,7 

III. EXPERIMENTAL VALIDATION 

Experimental verification of our calculated results is 
very important, because the geometrical optics assump- 
tions which underlie the ray-tracing methods are not well 
satisfied by the plasma plume. We are therefore planning to 
experimentally measure the effect of the plume on the re- 
flector antenna pattern. The issue for satellite communica- 
tions is the plume effect on the far-field radiation pattern. 
This is difficult to simulate in the laboratory, however, 
particularly with the arcjet vacuum requirements. A near- 
field experiment is more appropriate, and also provides a 
very sensitive test of our code. If the near-field calculations 
are verified correct, then our far-field predictions should be 
very accurate. 

The experimental configuration is shown in Fig. 2. A 
40-cm-diam reflector antenna directs collimated radiation 
at 10 GHz through the plasma plume. The diagnostic con- 
sists of a small movable rf probe coupled to a Hewlett- 


2978 Rev. Sci. Instrum. 61 (10), October 1990 0034-6748/90/102978-03502.00 © 1990 American Institute of Physics 2978 


67 


OF POOR QUALITY 


density 
is well 
o the 

ig ien- 

:ontour 
’-ast- 
rc ave 
cfor of 
er, and 
e< jm. 
ah^e of 
eled as 
ir re- 
st iral 

ite the 
le and 
n\ ted 

lh the 
lc >at- 
Uns- 
irmed, 
g <m- 

>i ing 

Deteri- 
er ob- 


u! is 
sump- 
■t well 
ii to 
h_re- 
jnica- 
.r *n. 
v< sr, 
near- 
ides a 
it ns 
iltr-be 

: A 
a_)n 
con- 
d-t- 



FIG. 1. Contour plot of the arejet electron density profile. Measurement 
locations are denoted by X . The contours are given by the least-squares 
fit n f — 9.1 X 10 M exp( — 0.0190)/^ cm -3 (0 in degrees, r in cm). 

Packard 8510 network analyzer. With the arejet operating 
the antenna radiation is deflected upward, since the index 
of refraction is less than unity within the plasma. This 
changes both the magnitude and phase of the near-field 
distribution, as shown in Fig. 3. These calculations are for 
the probe location shown as a dotted line in Fig. 2, with 
Y = 0 referenced to the arejet nozzle. The effect is rela- 
tively small, and concentrated toward the region closest to 
the plume. Here the maximum induced phase shift is —8° 
and the magnitude change the order of 2.5 decibels. 

The arejet thruster requires a modest vacuum for op- 
eration, typically — 10~ 4 Torr. Microwave reflections, par- 
ticularly from the wall behind the rf probe, are a consid- 
erable problem because of the relatively small tank 
dimensions ( — 1 mXl mXl m). Some attenuation will be 
obtained with vacuum compatible microwave absorbing 
material (Eccosorb FDS), which provides —3 dB/mm at- 
tenuation at 3 GHz. The absorbing material will damp 
multiple reflections. 

To eliminate primary reflections we will make use of 
the time domain gating feature of the HP8510 analyzer. In 
this technique measurements are made over a broad fre- 



FIG. 2. Schematic of the measurement configuration, showing approxi 
mate ray trajectories in the presence of the plasma. 
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FIG. 3. Calculated near-field magnitude and phase (referenced to the 
antenna plane) with and without the plasma. 

quency range in the ac steady state, and Fourier trans- 
formed to the time domain. Gating can then be applied in 
the time domain to eliminate the reflections, since they are 
shifted in time by the extra path length traveled. The gated 
data are then transformed back into the frequency domain. 
Because the measurements do not extend to zero frequency 
(co>co p ) t the FFT is performed over positive frequencies 
only. This operation preserves distinct time domain fea- 
tures (uniform response in the frequency domain trans- 
forms to a pulse in the time domain), although the time 
domain data are complex. 

Several sources of error must be considered when using 
time gating. 8 The finite sampling used for the FFT causes 
aliasing if the extra path length for reflections is too large. 
If the path difference is too small, the reflections will not be 
time separated from the main signal ( which iis broadened 
by dispersion and the limited range of frequency measure- 
ments). Aliasing is not expected to be a problem, because 
of the small tank dimensions and use of some absorbing 
material. Separation of the reflections may be difficult, 
however. 

The plasma is dispersive and distorts the transformed 
time domain pulse shape, complicating the task of time 
gating. A numerical simulation is shown in Fig. 4, where 
the frequency spectrum of the input signal is assumed to be 
Gaussian with a center frequency of 15 GHz and a FWHM 
of 2.94 GHz. The magnitude of the transformed pulse 
shape of the input Signal is shown in Fig. 4(a). The cor- 
responding FWHM of the pulse is 0.3 ns. The output of 
this signal after propagating through a homogeneous 
plasma slab of thickness d — 20 cm is computed, by assum- 
ing that each fr equency co mponent undergoes a phase shift 
of — ( cod/c ) yjl—col/o 1 . This is shown in Fig. 4(b) for 
different plasma frequencies. It is observed that as long as 
the frequency range of measurement is sufficiently above 
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(a) Input Pulse 



(b) Simulated Output Pulse Shape 



FIG. 4. Simulated time domain output of (a) the input pulse shape and 
(b) the output pulse shapes after propagation through a 20 cm plasma 
slab. 

the plasma frequency, no significant pulse distortion will 
occur. The observed pulse delay is also not significant. 

We have set up a test facility consisting of an antenna 
chamber, antennas, and HP8510 analyzer. We have veri- 
fied basic system operation, using time gating on induced 
reflections. We have successfully compared predicted and 
measured effects of various dielectrics (such as paraffin) to 
high accuracy. Our current facility does not have vacuum 
capabilities, however. After completing this proof-of- 
principle work, we plan to set up the measurement config- 


uration at NASA Lewis on their arcjet vacuum facility. 2 
This is expected to begin during the summer of 1990. 

IV. CONCLUSIONS 

Electric satellite positioning and propulsion is now 
practical because of the substantial power supplies now 
designed into the spacecraft. One of the most critical issues 
for its use on communication satellites is the effect of the 
thruster plasma on microwave communication systems. 
We have done extensive modeling and analysis of this 
problem, and are now verifying this work experimentally. 
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ABSTRACT 

The effect of an arcjet plume on the performance of satellite reflector antennas is 
studied. The arcjet plume is modeled as a weakly ionized plasma. Spatial permittivity 
distribution of the plume is approximated using the available measurement data of the 
electron density profile and a cold plasma model. Geometrical optics is applied to 
determine the ray paths as well as the transmitted fields through the inhomogeneous plume. 
The ray optics results are compared against several exact solutions for scattering from 
inhomogeneous dielectrics and good agreement is observed for sufficiently large scatterer 
size. The far field antenna patterns of the reflector in the presence of the plume are 
calculated from the transmitted ray fields using a ray-tube integration scheme. For arcjet 
prototypes in the 1-kW class, the plume effect on the antenna performance is small. As the 
electron density increases (corresponding to higher arcjet power), a gradual degradation of 
the main beam and sidelobe level is observed. In addition, the main beam tends to squint 
away from the plume region. 


This work is supported by the NASA Lewis Research Center under Cooperative 
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70 


1 INTRODUCTION 

An arcjet thruster is an electrothermal propulsion device being developed by NASA 
and the aerospace industry for stationkeeping and attitude control on future generation 
satellites. Because of its high efficiency, this device could potentially extend the lifetime of 
modem communication satellites up to ten years. The operating principle of the arcjet 
consists of heating a gas mixture (e.g., nitrogen-hydrogen) by means of a high temperature 
arc discharge. The heated propellant is then ejected through a nozzle to produce thrust. 
Recent arcjet research indicates that specific impulse levels well above 400 sec are attainable 
at the 1-kW level [1], making arcjets a very attractive candidate to replace the existing 
resistojet technology for auxiliary propulsion on satellites. One important issue which 
arises in the actual spacecraft integration of arcjet technology is the environmental impact of 
the arcjet plume, especially on microwave communication systems. The arcjet produces a 
weakly ionized plasma, which modifies the free space normally surrounding the satellite 
and causes refraction and attenuation of the microwave signals to and from the spacecraft. 
This is a primary user concern that must be addressed quantitatively before the arcjet 
technology can be fully accepted by the communication satellite industry. 

In this paper, the effect of an arcjet plume on the performance of on-board reflector 
antennas is studied (Fig. 1). The arcjet plume is modeled as an inhomogeneous scatterer 
with a smoothly varying index of refraction. The spatial distribution of the refractive index 
is related to the measured electron density data [2] based on the cold plasma model. To 
study the propagation of the microwave signal through the inhomogeneous plasma plume, 
geometrical optics is used. Choosing geometrical optics to attack this problem is motivated 
by several considerations: (a) the plume size of interest is more than several wavelengths in 
extent, (b) the permittivity distribution of the plume is smoothly varying, and (c) a ray 
approach is compatible with existing reflector analysis codes. A "shooting and bouncing 
ray” (SBR) scheme, developed earlier for characterizing the radar cross section of engine 
inlets [3], is adopted in implementing the present analysis. A set of rays is shot from the 
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reflector aperture toward the plume. The trajectory, phase, amplitude and polarization of 
the high frequency ray field transmitted through the plume are computed numerically based 
on the laws of geometrical optics [4]. To validate the geometrical optics results, 
comparisons with several exact solutions to scattering from inhomogeneous dielectric 
bodies are made. Once the ray optics field transmitted through the plume is found on a 
near-field aperture, the far field performance degradation of the reflector in the presence of 
the plume can be fuUy characterized. Several plume-to-reflector configurations are studied. 
The resulting gain loss, pointing error and sidelobe degradation are crucial in determining 
realistic constraints on the location and maximum power output of the arcjet for spacecraft 
integration. 


T PLUME MODEL 


Available experimental data taken from Langmuir probe measurements indicate that 
the arcjet plume consists of a weakly ionized plasma with electron density which vanes as a 
function of the operating conditions and the distance from the nozzle [2]. The measured 
electron temperature of the plasma is well below 1 eV, indicating that a cold plasma model 
should be adequate in describing the plume. It is also estimated that both the collisional 
effect and the anisotropy due to the earth magnetic field are very small [5]. Therefore, the 
refractive index for this unmagnetized cold plasma depends only on the electron density and 

can be modeled as 


N 2 


e_ 

eo 



2 n e e 2 
co = — 5 — 
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where c * p is the plasma frequency and n e is the electron density. From experimental data, 
the spatial distribution of the plume electron density is observed to obey approximately the 
functional form 
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where r is the distance from the nozzle and 0 is measured from the plume center line along 
the x-axis (see Fig. 2(a)). Numerically, a! is equal to the on-axis electron density 1 cm 
away from the arcjet nozzle. It is a strong function of the arcjet power, but also depends on 
other parameters such as propellant type and mass flow rate. Typical values for the 
constants aj and a are found, through curve-fitting to experimental data, to be 3xl0 12 cm" 1 
and 0.03 deg 1 , respectively. A more detailed description of the curve-fitting algorithm to 
the experimental data can be found in [5]. Using the above values, contour plots of the 
relative dielectric constant are generated. Figs. 2(b) and 2(c) are the plume profiles at 3 
GHz and 1 GHz, respectively. The contour levels are 0.02 apart with the outermost 
contour representing a relative dielectric constant of 0.98. Subsequent inner contours have 
decreasing dielectric constants. Very close to the nozzle, the dielectric constant can become 
negative due to the high electron density (and thus high plasma frequency) there. This 
corresponds to the cutoff for wave propagation (gkc o p ). Comparing the two contour plots, 
we observe that the physical size of the plume increases as the frequency of the microwave 
signal is decreased. On the other hand, the electrical size of the plume (measured in terms 
of the free space wavelength) remains essentially independent of the frequency. 


3 RAY TRACING THROUGH I NHOMOGENEOUS PLUME 
Once the plume model is established, geometrical optics is applied to calculate the 
propagation of the microwave signal through the plume. The solution of geometrical optics 
fields in an isotropic inhomogeneous medium was first derived by Luneburg more than 
thirty years ago and is documented in detail by Kline and Kay [4], Bom and Wolf [6] and 
most recendy, Bremmer and Lee [7]. Most numerical implementations of the ray solution 
for arbitrary media have been associated with studies of radio wave propagation in the 
ionosphere [8] and have not been widely reported in the antenna community. Senol and 
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Romine [9] used ray optics to study the effect of rocket plume on antenna performance, but 
their work was based on a scalar theory. Below the relevant equations for determining the 
ray trajectories as well as the associated phase, amplitude and polarization of the field along 
the ray are described. The detailed numerical implementation of these equations can be 
found in [10]. The trajectory of a ray in space is governed by the vector differential 
equation 

i(Nf) = VN (3) 

where t is a unit tangent vector of the ray at a point in space and the arclength s is measured 
along the curved ray trajectory (see Fig. 3). For an arbitrary inhomogeneous medium with 
index of refraction N(r*), the trajectory of a ray can be determined by solving the above 
equation numerically. For example, by finite differencing, we obtain the difference 
equation 

NC? +A? ) t (t + Ar*) = N (t) t(r) + As VNC?+ At/2 ) (4) 

which allows us to solve for the subsequent ray trajectory once the initial direction t of a 
ray at point t is given. When the medium is homogeneous, the ray trajectory reduces to a 
straight line, as expected. From the above expression, we also observe that the ray 
trajectories are continuously being altered in the direction of the refractive index gradient. 
Therefore, the ray paths will tend to bend toward regions of higher refractive index. Since 
the refractive index of the plasma is always less than unity, it is expected that the rays will 
bend away from the plume. 

Once the ray trajectory is determined, the phase variation along the ray path is given 
by Fermat's relationship: 

P2 

[S(Pi) - S(P 2 )] = J N ds (5) 

Pi 
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The complete geometrical optics field along the ray can also be determined. The field at 
point P2 on the ray path is related to that at Pi by the relationship 


fi(P 2 ) = e 2 E(Pi)VN] 7 N" 2 (DF) exp{-jko[S(P2)-S(Pi)] } (6) 

The quantity DF, or the so-called divergence factor, can be interpreted physically as a 
factor which accounts for the spreading of the differential ray tube such that the energy 
carried by the ray tube remains constant (Fig. 3 ). The divergence factor is derivable from 
the first order transport equation and is given by 

P2 

DF = exp { - i J (V • t ) ds } ( 7 a) 

Pi 

or 

P2 

DF = exp { - 5 J ( ^ + gj) ds ) ( 7 b) 


where R! and R 2 are the two principal radii of curvature of the wavefront along the ray. 
Because of its simple geometrical interpretation, equation ( 7 b) is better suited for numerical 
implementation of the divergence factor than ( 7 a). 

A A 

The polarization vector e 2 is related to the polarization vector at Pj, ej, in the 
following manner: 

e 2 = n cos(0q-0) + b sin(0Q-0) (8) 


where 
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In the above expressions, n , b and t are respectively the normal, binormal and tangent 
vectors along the ray (see Fig. 4). t is the "torsion" of the ray trajectory in space. It is a 
measure of the turning of the binormal vector. Note that the polarization vector depends on 
the torsion of the ray trajectory. For planar ray paths, the torsion is zero and therefore the 
electric field remains at a constant angle with respect to the normal and binormal vectors 
along the ray. 

Based on the ray optics equations, a two-dimensional and a three-dimensional ray- 
tracing program have been developed for arbitrary inhomogeneous media. Results of the 
programs have been verified against the exact solution in "Maxwell's fisheye" which has a 
refractive index described by N(r) = Ng / [ 1 + (r/a) 2 ]. The circular ray trajectories 
generated by the numerical routine agree very well with the exact solution. The divergence 
factor calculation is also checked against the exact solution derived by Bremmer and Lee [7] 
and found to be accurate. Shown in Fig. 5 is ray trajectory generated using the three- 
dimensional ray code for a "two-dimensional fisheye" medium which has refractive index 
profile N(p) = N 0 / [ 1 + (p/a) 2 ]. With a non-zero z-component in the initial ray direction, 

the ray trajectory spirals about the z-axis in a helix. Also shown as tick marks are the 
electric field vectors along the ray path. 

4 VALIDATION OF RAY RESULTS 

In order to validate the geometrical optics theory and its numerical implementation, 
several scattering problems are investigated. First, the two-dimensional scattering of a 
plane wave by an inhomogeneous circular cylinder with a continuous and radially-varying 
refractive index is considered (see Fig. 6). To establish a reference solution, the series 
expansion solution to scattering from a multi-layered circular cylinder [11] is implemented 
to approximate the exact solution. A smoothly varying refractive index profile 
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N(p) = WP) = 


p < a 


(9) 


*^b [1 - c cos ( “ ')] » 


l 1 , 


p > a 


is chosen, where b=H + e r (0)]/2 and c=[l- £,(())]/[ l+e r (0)]. To generate a numerical 
solution utilizing the two-dimensional ray tracing routine, a set of collimated rays from the 
incident plane wave direction are traced through the inhomogeneous dielectric. Shown in 
Fig. 7 is the ray trajectories inside the inhomogeneous cylinder for e/0^0.5. The rays 
tend to bend away from the center region which has a lower refractive index. Once the ray 
fields entering and exiting the scatterer are determined, the scattered far field can be found 
by carrying out the radiation integral one ray tube at a time. This "ray-tube integration 

formulation is similar to those reported earlier [3],[12]. 

Figs. 8(a)-(d) are comparisons of the E-polarized scattering widths of a cylinder with 
ej(0)=0.5 calculated using the ray formulation and those obtained from the multi-layered 
solution for cylinders of different radii. Ten layers per wavelength is used in the multi- 
layered solution in order to simulate a smooth refractive index profile. It is observed that 
for scatterer diameter above 3 free space wavelengths, the agreement between the two 
results is very good in the forward scattering region. This is as expected since the 
geometrical optics theory should become more accurate as the frequency is increased. For 
angles far away from the forward(4>=0°) direction, the exact solution is more than 30 dB 

below the forward scattering strength. In this region, the accuracy of the ray results is 
limited by the numerical noise floor arising from the finite ray density used in the 
calculation. The same comparisons can be made for the H-polarized case. In fact, based 
on the ray formulation, it can be shown that the scattering pattern of the H-polarization is 
identical to that of the E-polarized case. From the exact solution, it is also found that as the 
scatterer size is increased beyond 3 X, the scattering widths of the two polarizations become 
essentially identical. This is consistent with the commonly held notion that the scattering 
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pattern is not very sensitive to polarization in the high frequency regime. Fig. 9 shows the 
comparison of the scattering widths of a cylinder of diameter 6 A calculated using the ray 
approach and those obtained from the multi-layered solution for e r (0)=3.0. In this case the 
rays are attracted toward the center and caustic points exist inside the cylinder. Again the 
ray results agree very well with the multi-layered solution. 

Next, the scattering from an inhomogeneous sphere is earned out as a three- 
dimensional validation example. A refractive index profile similar to (9) is chosen, except 
that the variable p is replaced by the three-dimensional radial varible r in the expression. 
The series expansion solution to scattering from a multi-layered sphere [13] is implemented 
to approximate the exact solution to the problem. The ray optics results are generated from 
the three-dimensional ray tracing code and a three-dimensional ray-tube integration scheme 
similar to that used in the two-dimensional case. The comparisons between the two results 
for a sphere of diameter 6 A are shown in Figs. 10(a) and 10(b) for £j.( 0)=0.5 and 3.0, 
respectively. Agreement between the ray results and the multi-layered solution is consistent 
with the two-dimensional cases. 

5. REFLECTOR PERFORMANCE DEGRADATION 

With the ray optics method validated, the effects of the arejet plume on reflector 
performance are investigated. (Some preliminary results have been reported earlier in [14] 
and [15]). Shown in Fig. 1 1(a) is the three-dimensional geometry of the case under study. 
The reflector is assumed to have a diameter of 60 cm (20 A at 10 GHz) and a 10-dB 
aperture edge taper. The initial ray bundle after reflection from the reflector surface is 
assumed to be perfectly collimated and has perfect linear polarization in the y-direction. 
Also shown in the figure are ray trajectories through a plume with parameter a^lxlO 14 
cm 1 . As expected the rays tend to bend away from the high electron density region near 
the nozzle. The direction, phase, amplitude and polarization of the ray field transmitted 
through the plume are computed on an exit aperture plane Za=50 cm. The far field of the 
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reflector in the presence of the plume can be related to the transmitted field on the aperture 
plane through the standard radiation integral as follows: 


— » e -jkr A A 

E = ^r- (0A e + 


( 10 ) 


A e 


L a* 


-A J I exp[ jk (cos<|) sin© x + sin<j> sin© y + cos© Zj,)] dxdy 
4 k 


z=z a plane 


- a 


A 

L 0 


x E ( r ) + T| 




L* 


— ► — > , A 

xH (r) } • z 


To evaluate the above aperture integral, we follow a "ray-tube integration" procedure 
similar to those used in Sec. 4 and sum up the contributions from all rays arriving on the 
aperture [3], [12]. Carrying out the ray-tube integration over the aperture plane results in 

A e * = ( ) X (Area) pro j exp[ jk(cos<J> sin© x* + sin<> sin© yj + cos© Za)] 

4k ith ray 

• Be, * • S(0, <|>) (ID 

where 

Be = [ sin<J> E y i + cosf E^j ] + [ cos© cos<J> H y>i - cos© sin<> H x>i ] 

= [ cos© cos<J) E y>i - cos© sin<J> E^j ] + r\ [ - sin<|> H y>i - cos<}» H x>i ] 
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= the projected area of the ray tube on the aperture plane, 
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(Area)in is the initial ray tube size before the plasma. 

In deriving (11), we have assumed that the ray direction is given byt = t x ^ + t y 9 + t z z, 
and that the magnetic field is related to the electric field through H = (txE )/Ti. 

The reflector patterns along the two principal planes for the geometry in Fig. 1 1 (a) are 
computed and plotted in Figs. 1 1 (b) and (c). Small but noticeable degradations in the main 
beam and the sidelobe regions are observed. In particular, the main beam is found to 
squint away from boresight in the H-plane. This effect can be interpreted physically if we 
ignore, to the first order, the refraction of the rays through the plume. We can then argue 
that the rays traveling closer to the nozzle must experience a smaller total phase shift than 
those farther away from the nozzle due to the lower permittivity of the inner plume region. 
Consequently, the transmitted field on the exit aperture z^ has a tilted phase distribution 
which leads to the main beam squinting away from the plume region. In this example, the 
beam (beamwidth=3.0°) experiences a 0.75° (0.25 BW) beam squint away from the 
boresight. This translates into a 0.94 dB gain loss. For very narrow beam antennas, this 
effect could become more serious. The E-plane pattern in Fig. 11(c) shows a higher 
sidelobe level than the H-plane pattern, and a noticeable cross polarization level due to the 
plume. Figs. 12(a) and (b) show the reflector far-field patterns under a higher power arcjet 
operating condtion ( ai =3xl0 14 cm' 1 ). The deterioration of the reflector performance is 
more severe. The beam squint and boresight gain loss are respectively, 1.5 and 4.84 dB. 
Finally, it is worth noting that the values of 3l used in the above examples are significantly 
higher than the measurements reported in [2], These larger values, through extrapolation 
of the experimental data, should correspond to arcjets operating in the 3-kW range. For 
ai =3xl0 12 cm 1 , corresponding to the 1-kW class laboratory prototypes used in [2], the 

reflector pattern degradation is very small. 

Next, a different plume configuration is studied where the plume is oriented 
horizontally along the y-axis and the nozzle position is located at (x=0, y=-30 cm, z=0), as 
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shown in Fig. 13(a). The reflector geometry and the arcjet power are kept identical to those 
of Fig. 12. The reflector patterns in the two principal planes are plotted in Figs. 13(b) and 
(c). In this case, the plume distribution has no symmetry along either the H- or the E-plane 
of the reflector and the beam squint phenomenon can be observed in both planes. In the H- 
plane, the beam squints toward the +6 direction, as expected. In the E-plane, it may first 
appear that the beam should squint away from the nozzle and toward the +0 direction. 
However, closer examination of the plume profile reveals that there is more plume blockage 
in the y>0 half of the aperture. Therefore, the beam squints toward the -0 direction. The 
gain loss at the boresight is 1.18 dB. By comparing the antenna patterns in Fig. 12 to 
those in Fig. 13, we notice that the overall pattern degradation is much less severe in the 
case presented in Fig. 13. We therefore conclude that the plume effect on reflector 
performance is strongly dependent on the actual plume-to-reflector configuration. 

6. SUMMARY 

The effect of an arcjet plume on the performance of satellite reflector antennas has 
been investigated. The arcjet plume is modeled as a weakly ionized plasma. Spatial 
permittivity distribution of the plume is approximated based on the cold plasma model and 
the available measurement data of the electron density profile. Geometrical optics is applied 
to determine the ray path as well as the transmitted field strength through the 
inhomogeneous plume. The scattered fields calculated from geometrical optics have been 
verified against the multi-layered solutions to radially inhomogeneous two-dimensional 
cylinders and three-dimensional spheres. The far field patterns of the reflector in the 
presence of the plume have been computed for various cases and the resulting degradations 
have been studied. As the electron density is increased (corresponding to higher arcjet 
power), a gradual degradation of the reflector pattern is observed. The most significant and 
interesting effect is a beam squinting phenomenon which results in antenna gain loss and 
pointing error. The amount of degradation is also shown to be strongly dependent on the 


81 


plume- to-reflector configuration. The capabilites developed in this work can be applied to 
model more realistic antenna configurations on satellites. The quantitative information 
generated will be crucial in determining realistic constraints on the location and maximum 
power output of the arcjet for spacecraft integration. 




The authors would like to thank L. M. Carney for supplying the electron density 
measurements used in this work and for helpful discussions on the plume characteristics. 
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(a) Coordinate geometry for describing the plume, (b) Relative dielectncconsfcmt 
profile of the plume at 3 GHz, trnd <c) 1 GHz. The contour levels are 0.02 apart. 





|DF | 2 


Differential Ray Tube Area at 


Differential Ray Tube Area at P 2 


Fig. 3. Curved ray trajectory in an inhomogeneous medium. Ray tube spreading is 
accounted for by the divergence factor (DF). 



Fig. 4. 


Three unit vectors n (normal), b (binomial) and t (tangent) can be defined along 
the ray trajectory in space. 
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Fig. 6. Two-dimensional scattering of a plane wave by an inhomogeneous circular 
cylinder, (a) Problem geometry, (b) Refractive index profile. 
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(a) 



(b) 



Fig. 8. Comparison of the scattering width of an inhomogeneous cylinder with 

calculated using ray optics and that obtained from the multi-layered 
solution, (a) 2a = IX. (b) 2a = 3X. 
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(C) 



(d) 



Fig. 8(Cont'd). (c) 2a = 5 \. (d) 2a = 7X. 
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Fig. 9. Comparison of the scattering width of an inhomogeneous cylinder with 
2a= 6X and e r (0)=3.0 calculated using ray optics and that obtained from 
the multi-layered solution. 
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Fig. 10. Comparion of the scattering cross section of an inhomogeneous sphere of 

diameter 6X calculated using ray optics and that obtained from the multi-layered 
solution, (a) Er(0)M).5. (b) £ r (0)=3.0. 
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(a) 



Fig. 11. Reflector pattern degradation due to a plasma plume (apslxlO 14 cm* 1 ). 

(a) Reflector and plume configuration. The reflector has a diameter of 60 cm 
(20X at 10 GHz) and an assumed aperture taper of 10 dB. The polarization is 
linear in y. (b) H-plane pattern, (c) E-plane pattern. 
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Fig. 12. Reflector pattern degradation due to a higher power arcjet plume (a^xlO^cnr 1 ). 

(a) H-plane pattern, (b) E-plane pattern. 
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(a) 



10 GHz 

a 1 =3xl0 14 cm- 1 


Fig. 1 3. Reflector pattern degradation due to a different plasma plume configuration 

with aj=3xl0 14 cm' 1 , (a) Reflector and plume geometry. The reflector has a 

diameter of 60 cm (20X at 10 GHz) and an assumed aperture taper of 10 dB. 
The polarization is linear in y. The plume is oriented horizontally along y and 
the nozzle position is located at (0, -30cm, 0). (b) H-plane pattern, (c) E-plane 
pattern. 
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HIGH FREQUENCY ELECTROMAGNETIC SCATTERING FROM AN 
INHOMOGENEOUS DIELECTRIC BODY BY RAY TRACING 


Hyeongdong Kim and Hao Ling 


Department of Electrical and Computer Engineering 
The University of Texas at Austin 
Austin, TX 78712-1084 


Abstract - A ray formulation is presented for treating the electromagnetic scattering 
from electrically large, inhomogeneous objects. A dense grid of rays representing the 
incident plane wave is shot toward the inhomogeneous object. At the scatterer boundary, 
reflected rays and refracted rays are generated due to discontinuity of the medium 
parameters. The trajectory, amplitude, phase and polarization of the refracted rays are 
traced inside the inhomogeneous object based on geometrical optics. Whenever the rays 
cross the scatterer surface, additional reflected/reffacted rays are generated and are tracked. 
This process is repeated until the intensities of the refracted/reflected rays become 
negligible. The contributions of the exiting rays to the total scattered field arc calculated by 
using the equivalence principle in conjuction with a ray-tube integration scheme. The ray 
formulation is applied to calculate the bistatic scattering from cylinders and spheres. The 
numerical results of the ray formulation are compared with exact series solutions for 
homogeneous and inhomogeneous objects. 


This work is supported by the National Science Foundation under Grant ECS-8657524 and 
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I. INTRODUCTION 


The study of electromagnetic scattering from inhomogeneous objects has received a 
great deal of attention because of its importance in many practical areas such as dielectric 
lens design, radar signatures from aircraft and missile plumes and power deposition in 
biological tissues. Several numerical techniques are available to treat such problems. 
Among them are the method of moments solution using either surface integral [1] or 
volume integral equations [2], the unimoment method [3], and hybrid methods [4], [5] 
combining the finite element method with the integral equation approach. All of these 
"numerically exact" methods, however, are limited to the low frequency regime where the 
electrical size of the scatterer is small, due to the computation time involved and the amount 
of computer memory required. For electrically large scatterers, approximate methods such 
as high frequency ray techniques become more attractive. 

Kouyoumjian et al. [6] used a modified geometrical optics (GO) approach to study 
the backscattering from electrically large, homogeneous dielectric objects. By tracing rays 
analytically inside homogeneous objects of special shape, good results were obtained in the 
backscattering direction. Chen [7] formulated the problem of bistatic scattering from 
homogeneous objects by considering diffracted rays as well as GO rays. However, no 
numerical results were reported in his work. In this paper, a ray approach is presented to 
characterize the scattering from electrically large, inhomogeneous objects. This work is 
motivated by our earlier study of space plasma effects on microwave reflector antenna 
performance [8]. Since the inhomogeneous plasma plume from an arcjet thruster is large in 
terms of the wavelength of the reflector antenna signal, ray theory is applicable. The ray 
tracing technique documented in this paper is a natural extension of the so-called "shooting 
and bouncing ray" (SBR) scheme, developed earlier for calculating the radar cross section 
of partially open inlet structures [9], A dense grid of geometrical optic rays representing an 
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incident plane wave is shot toward the inhomogeneous object. At the scatterer boundary , 
reflected rays and refracted rays are generated due to discontinuity of the medium 
parameters (see Fig. 1). The trajectory, amplitude, phase and polarization of the refracted 
ray are traced inside the inhomogeneous object based on the laws of geometrical optics. 
Whenever the rays cross the scatterer surface, additional reflected/ refracted rays are 
generated and are tracked. This process is repeated until the intensities of the 
refracted/reflected rays become negligible and the ray tracing algorithm can be terminated. 
The contributions of all the exiting rays to the exterior scattered field are then calculated by 
using the equivalence principle in conjuction with a ray-tube integration scheme [10]. 

The organization of this paper is as follows. In Sec. 2, the ray formulation is 
presented. First, the propagation of the GO field in an arbitrary inhomogeneous medium 
and the GO reflection/refraction problem at a curved dielectric interface are discussed. 
Once information on the exiting rays are found, a ray tube integration scheme, which 
consists of summing the scattered fields from individual ray tubes, is considered in order to 
arrive at the scattered field in the exterior region. In Sec. 3, numerical results for various 
inhomogeneous objects are presented. They are compared against the available exact 
solutions. Finally, some conclusions are summarized in Sec. 4. 
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II. FORMULATION 


The geometry of the scattering problem is shown in Fig. 1. An inhomogeneous 
object, which is defined by the constitutive parameters e and |i (or refractive index 
N=(e r ji r ) 1/2 ), is illuminated by an incident plane wave in free space. The refractive index of 
the scatterer is assumed to be smoothly varying, except at the scatterer boundary. The 
objective is to determine the scattered field using a ray tracing approach. To represent the 
incident plane wave, a large number of rays are shot toward the inhomogeneous object. The 
incident rays encounter the scatterer boundary and reflected rays and refracted rays are 
generated, e.g., at point 1 in Fig. 1. The refracted rays propagate inside the object and 
eventually encounter the scatterer boundary (point 2 in Fig. 1). Again reflected rays and 
refracted rays are generated, except this time the refracted rays exit the scatterer while the 
reflected rays bounce off the boundary and propagate inside the inhomogeneous object. The 
intensities of the reflected rays should become weaker after every bounce against the 
boundary and this process can be terminated after several bounces. To account for the ray 
reflection and refraction at the curved scatterer boundary, the geometrical optics solution 
described in [4] is summarized in Sec. 2.1. Next, ray tracing in an arbitrary 
inhomogeneous medium is implemented numerically. The relevant equations [3] for the ray 
trajectory, the associated phase, amplitude and polarization of the ray field are detailed in 
Sec. 2.2. Finally, whenever a ray exits the scatterer boundary, the contribution of the ray to 
the total scattered field is calculated by a ray- tube integration scheme [9,1 1], The integration 
is taken over the projection of the ray tube on the scatterer boundary and the fields on the 
scatterer boundary is approximated by the ray fields. This scheme is described in Sec. 2.3. 
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2.1. GO Reflection/Refraction at the Curved Scatterer Boundary 

When a ray is incident upon the scatterer boundary, a reflected ray and a refracted 
ray are generated due to the discontinuity in the medium parameters. To determine the 
fields associated with the reflected and refracted rays as they propagate away from the 
boundary, their initial values including the ray directions, amplitudes, polarizations and 
wavefront curvatures are needed at the impact point. Let tj, t r , tt b® the direction vectors 

of the incident, reflected and refracted rays, respectively. The plane defined by tj and n is 

A A 

the plane of incidence and v is a unit vector perpendicular to that plane (see Fig. 2). t r> tt 
can be determined by Snell's law as follows: 

tr • U ) = tj • n , t r • V = 0 , 1 tr I = 1 

m ( t r • n ) = n 2 ( ti • n ) , t t • v = 0 , | t t | = 1 ^ 


For the amplitude and polarization calculation, Fresnel's formulas for a planar interface can 
be used for a curved interface under the high frequency approximation. The incident field 
is first decomposed into its perpendicular and parallel components as follows (refer to Fig. 


2 ): 


E 1 = ( E'*<t> )<))+( E *0 )6 


( 2 ) 


The reflected and transmitted fields at the impact point on the boundary are then obtained by 
using the planar reflection and transmission coefficients: 

E r = R( E^O? + R( E 1 *©)? (3) 

• 

E 1 = T( E 1 ^ )<j> + T( E 1 *© )0 

where 
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^1 ^ ^ 

0 = ucos 0i - nsin 0, 

'T 

0 =-ucos 0i - nsin 0i 

0 =-ucos 0 t - nsin 0 t 

-s ^ 

<J) = 0 = <{> = V 

01 and 0 t : incident angle and transmitted angle. 

Information on the reflected and refracted wavefronts are necessary in order to 
determine the ray fields away from the scatterer boundary. The curvatures of the reflected 
and refracted wavefronts can be compactly represented by the wavefront curvature 
matrices, Q" and Q* [16], from which the principal radii of curvature and the principal 
directions can be extracted. The curvature matrices Qf and Q* can be determined by 
enforcing the phase matching condition on the scatterer boundary. This results in the 
following matrix equations, relating with the incident wavefront curvature matrix Q and the 
boundary surface curvature matrix Q s , to and Q T : 

B?Q r B r = BjQ'B, -2cos 0iQ s (4) 
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= BiQ'Bi + ( ^cos 0 t -cos 0i)Q s 


where 


B ilt = 


0 

t 

/N 

•u 

« i,l ~ 
0 »v 


cos 0i, t 

0 

^i, 

<t> 

,t 

•u 

1 


0 

1 


0 *U 

<> 

1 


-COS 0; 

0 

•~r 

<(> *U 

<> 

1 


L 0 

1 . 


B r = 


The curvature matrices Q s , Q\ Q r , Q‘ are described by the basis vector sets (u, v ), 
( 0 1 , 'ty ) (0,<t>) ( 0 , <(> ), respectively. Once and Q* are determined from the above 
matrix equations, they may be diagonalized in the standard fashion to arrive at their 
eigenvectors (principal vectors) and eigenvalues (principal radii of curvature). These 
information are used subsequently to calculate the ray field inside the inhomogeneeous 


scatterer. 


2.2. Ray Propagation in Inhomogeneous Media 

The propagation of a ray field inside the inhomogeneous scatterer, including the ray 
trajectories and the associated phase, amplitude and polarization, is considered. The 
solution of geometrical optic fields in inhomogeneous media was first derived by Luneburg 
more than 30 years ago and has been documented in detail by Kline and Kay [13], Bom 
and Wolf [14] and most recently, Bremmer and Lee [15]. It is important to point out that 
the ray optics field is a high-frequency approximation to the exact full-wave solution. In 
referring to high frequencies, the key is the relation of the scatterer size and curvature of the 
wavefront to wavelength. Therefore, the following summarized results concerning ray 
tracing are valid only when the medium properties, represented by e and \i or N=(e^i r ) 1/2 , 
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vary slowly and continuously over distances of several wavelengths. When the variation 
of the medium is not smooth, a full wave analysis is required. Ray optics also fails near 
focal points, or on caustic surfaces, where it predicts an infinite field. 

The ray trajectory in space is governed by the vector differential equation: 


-d- ( nt ) = v n 

ds 


(5) 


where t is a unit tangent vector at a point along the ray. The arclength s is measured along 
the ray trajectory, and N is the refractive index. For an arbitrary inhomogeneous medium 
with a continuous index of refraction N, the trajectory of each ray can be traced by solving 
the above vector differential equation. At points where N is discontinuous, the differential 
equation does not apply. Instead, a reflected ray and a refracted ray are generated, as 
discussed earlier. 

Once the ray trajectory is determined, the phase variation along the the ray is given 
by Fermat's relationship: 

S(Pi ) = S(P 2 ) + f n ds (6) 

J pi 


The complete ray optics field at point P 2 can be determined from the values at point Pi by 
the following equation: 


E(P 2 ) = e 2 E(Pi) fnTTn^ (DF) exp{-jko[S(P 2 )-S(Pi)]} ( 7 ) 

The Field amplitude for a lossless medium is governed by the quantity (DF), the so-called 
divergence factor, which accounts for the expansion of the differential ray tube so that the 
energy carried by the ray tube remains constant (see Fig. 3). The divergence factor is given 
by 
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( 8 ) 


f 


DF = exp{-i- (^- + ^")ds} 


R 2 


where R x and R 2 are the two principal radii of curvature of the wavefront along the ray. 
The directions in which the radii of curvature at a point on the wavefront have minimum 
and maximum values define the two principal directions. 

The polarization vector e 2 at P 2 is related to the polarization vector ei at Pj by the 

equation: 

e 2 = n 2 cos(0 o -0) + b 2 sin(0 o -0) ( 9 ) 


where 


0 O = tan _1 (— ) , 0 = 
«o 



xds 


- - un (tXt'K' 

ei = na 0 + b(J 0 , t ■ — 

M 

In the above equation, n, b and t are respectively the normal, binormal and tangent vectors 
along the ray, and the prime represents the derivative with respect to the arclength s. x is 
the torsion of the ray paths in space. For planar ray paths, the torsion is zero and the 
electric field remains at a constant angle with respect to the normal and binormal vectors 
along the ray. 


2.3. Ray-Tube Integration 

Once the total fields on the scatterer boundary (E, H) are calculated, the scattered 
field in the exterior region can be calculated based on the equivalence principle [12]. As 
shown in Fig. 4, the equivalent surface currents are related to the surface fields by 
J s = n x H and M s = -n x Ewhere n is the outward surface normal vector. They radiate in 
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free space to give rise to the scattered field. The scattered far field can be calculated by the 
following radiation integral: 



where <j>, 0 are unit vectors in the direction <p, 0 and r’ is the position vector of the 
scatterer surface. However, the fields on the scatterer boundary (E, H), calculated by ray 
tracing, are not available on a uniform grid of points over the surface. Direct integration of 
the above radiation integral, therefore, cannot be easily carried out and the ray-tube 
integration scheme described in [9,11] is used. Since the incident plane wave can be 
considered as a superposition of small ray tubes, the contribution of each ray tube to the 
scattered field is calculated separately. Whenever the ray tube is reflected/refracted at the 
scatterer boundary, the ray tube contribution to the scattered field is calculated. 

Referring to Fig. 5, when a small tube of ray exits the scatterer, the field at ? on the 
scatterer surface can be approximated as: 

E(r) = E(r a ) exp { -jkt»(r-r a ) } (11) 

A — ► ^ 

H(r) = H(r a ) exp { -jkt»(r-r a ) } = D< ^— exp ( -jkt»(r-r a ) } 

A) 

A 

where ? a describes the position perpendicular to the ray direction and t is the direction of 
the exit ray. In other words, the field on the scatterer boundary is assumed to have the 
same magnitude as the ray tube and linear phase variation across the ray tube projection on 
the surface. The radiation integral in (10) can be applied to calculate the scattered field 
contribution of this small ray tube: 
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471 I 
1 J V 


ray tube 
projection 
on S 


ejk*r 


/r 

-4> 




e 


xE s (?a) + zd 


e 

<fr 


xHs(ra)|e-j kt * (r -ra)»n ds' 


(12) 


4k 


e 


xE s (r a ) + Zo 

/N 

e 

xH s (ra) • 


_ <t>. 

J 


rn (AA) S(0,<}>) 


where 


S ( 6, <f> ) = 


l 


ejk(k-t)-^ -r a ) ds’ 


(AA) J ra y ^be 

projection 


XN * 

A A = Area of ray tube projection on surface S = (Ray tube area / cos ( n, t )) 


This process is repeated as the ray tube undergoes multiple reflections with the scatterer 
boundary, until the remaining energy of the ray tube is negligible with respect to the energy 
carried by the incident ray tube. Finally, 



as the total scattered field is the superposition of the scattered field of all individual ray 
tubes 
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HI. NUMERICAL RESULTS 


Based on the ray formulation in Sec. 2, both a two-dimensional and a three- 
dimensional ray tracing code have been implemented. It is found from numerical 
experimentation that launching incident rays from a nonuniform grid (refer to Fig. 6) is 
more efficient computationally than using uniform ray density. A higher ray density should 
be used on the outer region of the scatterer so that the number of rays per unit perimeter (or 
in the three-dimensional case, per unit surface area) launched toward the scatterer boundary 
is fixed. Convergence test is performed and it is observed that more than four rays per 
wavelength of the scatterer perimeter are necessary for convergent results in the two- 
dimensional case. A representative convergence test for the backscattering from a circular 
cylinder is shown in Fig. 7. 

Since the ray formulation is inherently a high-frequency approximation, it is 
instructive to investigate the lower bound of validity of the ray results. The ray formulation 
is applied to the scattering from homogeneous circular cylinders of different sizes and the 
results shown in Figs. 8-10. Both the TM and TE bistatic scattering patterns generated 
using the ray formulation are superimposed with the exact series solutions for cylinders 
with diameter IX, 4X and &X. Note that the ray results, as expected, improve in accuracy 
as the size of the object increases. The lower limit of the scatterer size, to which the ray 
scattering formulation can be applied, is about three free-space wavelengths. It is also 
observed that while the agreement with the exact solution is very good for the TM 
polarization, the agreement is not as good for the TE case, particularly near 0 = 150°. It is 
not clear what causes the poor agreement, although the same observation can be made from 
Kouyoumjian's work [6]. 

The ray scattering formulation is next applied to two-dimensional radially 
inhomogeneous circular cylinders for both TM and TE incident waves. Since the exact 
series solution is available for radially multi-layered cylinders [18], the solution of the 
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radially inhomogeneous cylinder is obtained by modeling it as a cylinder of many layers 
with different refractive indices. More than 10 layers per wavelength are used to 
approximate the exact solution for the inhomogeneous cylinder. The ray formulation is 
applied to cylinders with a smoothly varying refractive index profile of the form: 



(13) 


where 

, 1+Er(0) , l-£r(0) 

2 ’ 1 + 6 ,( 0 ) 

Shown in Figs. 1 1 - 13 are the comparisons of the ray results with the layered solutions for 
different e^O) with d=4. The diameter of the cylinders is chosen to be 6 X.. It is observed 
that the ray results agree well with the exact solutions for the TM cases. For the TE cases, 
the agreements are not as close as those for the TM cases. This could be caused by the GO 
reflection/refraction modeling at the scatterer boundary also noted earlier in the 
homogeneous case. As a more general application, the ray formulation is applied to an 
inhomogeneous cylinder with angular as well as radial variations in the refractive index. 
Fig. 14 shows the bistatic scattering patterns. 

For three-dimensional cases, the ray formulation is applied to scattering from 
inhomogeneous spheres. The three-dimensional refractive index profile is assumed to have 
the same radial variation as that of the above two-dimensional inhomogeneous cylinder, 
except that the variable p in expression (13) is replaced by the three-dimensional radial 
variable r. For comparison, the exact solution is obtained by using the multi-layered 
model [19] of radially inhomogeneous spheres, as in the two-dimensional case. Shown in 
Fig. 15 are the E-plane bistatic scattering patterns for several three-dimensional scatterers. 
From the numerical results, it was found that the E-plane and the H-plane scattering 
patterns are almost the same, as expected for large scatterers. Fig. 15(a) shows a sphere 
without any discontinuity in the refractive index at the boundary. In this case, no reflection 


112 


is generated at the scatterer boundary and good agreement is obtained between the ray 
results and the exact solution. Figs. 15(b) and 15(c) are the results for, respectively, a 
homogeneous sphere and an inhomogeneous sphere with a sharp discontinuity at the 
boundary. The agreement between the ray results and the exact solutions is only good in 
the forward scattering region. It is also found that the agreement deteriorates as the 
refractive index of the scatterer increases. A possible reason for the discrepancies can be 
traced back to the reflection/refraction problem for the parallel (TE) polarization case. Since 
the three-dimensional results involve both polarizations, the same problem is manifested in 
the results. In addition, more accurate high-frequency analysis such as the geometrical 
theory of diffraction for materials should be included in the ray tracing. However, we can 
see that the approximate scattering patterns can be calculated by this ray formulation, and it 
is expected that the ray formulation will be more accurate for the larger scatterers. 

As a final note, the present ray approach uses GO fields on the scatterer surface. If 
caustic points are located on or near the surface, GO fields are not accurate and the scattered 
field calculated based on the surface GO fields will not be accurate. This problem can be 
solved by choosing a Huygens' surface larger than the scatterer boundary so that the 
caustic points can be avoided. Another problem in the present formulation is the 
occurance of total internal reflection, when a ray is incident at an angle greater than the 
critical angle in the denser medium. Total internal reflection at a curved interface is a 
complicated phenomenon and the conventional Fresnel formulas should be modified. More 
detailed discussions are found in [20] - [22]. 
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IV. CONCLUSION 


A ray tracing scheme was presented to calculate the bistatic scattering patterns for 
electrically large, inhomogeneous objects. Inside the inhomogeneous object, the reflected 
and refracted ray bundles from the scatterer boundary were traced based on the laws of 
geometrical optics. From the fields and the "footprints" of the ray bundles on the scatterer 
boundary, the scattered field pattern was obtained by a ray-tube integration technique. To 
test the validity and accuracy of the ray tracing scheme, it was applied to the 
electromagnetic scattering of both two-dimensional and three-dimensional objects. It was 
shown that the ray results agree approximately with the exact series solutions. In the 
forward scattering region, the agreement is very good. As expected, it was observed that 
the accuracy of the ray scheme improved as the size of the object was increased. From 
numerical experiments, it was found that the ray tracing scheme is valid for objects greater 
than about three free-space wavelengths. Also, it was observed that the accuracy of the ray 
scheme deteriorates for objects of high refractive index. More accurate account of the 
scattering mechanisms at the object boundary using the geometrical theory of diffraction 
may alleviate the limitations of this geometrical optics scheme. 
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Figure Captions 


Fig. 1. The geometry of ray scheme. 

Fig. 2. Local coordinate system of ray tube at the surface scatterer. 

Fig. 3. Calculation of divergence factor for amplitude variation. 

Fig. 4. The geometry of scattering problem, (a) original problem, (b) equivalent problem. 

Fig. 5.The calculation of the total field on the Huygen surface (at r) from the ray tube field 
(at rD - 

Fig. 6. Non-uniform ray density for effective convergence. 

Fig. 7. Convergence test of back scattering field as the number of rays per perimeter length 
in wavelength increases* 

Fig. 8. The scattering pattern of homogeneous circular with 1 X diameter, (a) TM. (b) TE. 

Fig. 9. The scattering pattern of homogeneous circular with 4 X diameter, (a) TM. (b) TE. 

Fig. 10. The scattering pattern of homogeneous circular with 8 X diameter, (a) TM. (b) 
TE. 

Fig. 11. The scattering pattern of inhomogeneous circular cylinder with 2a=6X, £r(0)=3 
and £r(a)=4. (a) TM. (b) TE. 

Fig. 12. The scattering pattern of inhomogeneous circular cylinder with 2a=6X, e r (0)=6 
and £r(a)=4. (a) TM. (b) TE. 

Fig. 13. The scattering pattern of inhomogeneous circular cylinder with 2a=6X,, £r(0)=8 
and er(a)=4. (a) TM. (b) TE. 

Fig. 14. The scattering pattern of angularly inhomogeneous circular cylinder. (a) TM. 
(b) TE. 

Fig 15. The scattering pattern of spheres, (a) inhomogeneous with discontinuous 
boundary, (a) homogeneous, (c) inhomogeneous discontinuous boundary. 
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Inhomogeneous 

Scatterer 


Fig. 1 The geometry of ray scheme 
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Differential Ray Tube Area at Pj 


Differential Ray Tube Area at P 2 


Fig. 3. Calculation of divergence factor for amplitude variation 
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Fig. 4. The geometry of scattering problem 

(a) Original problem, (b) Equivalent problem. 
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Huygens’ Surface ( S ) 

O 4 Ray Tube 

Fig. 5. The calculation of the field on the Huygens' surface 
from the ray tube field 



Nonuniform 
ray density 



A L< 0.1 X 


Fig. 6. Nonuniform ray density for effective convergence 



NO. OF RAYS PER PERIMETER 


Fig. 7. Convergence test of back scattering field of cylinders as the number 
of rays per perimeter length in wavelegth increases 
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Fig. 14. The scattering pattern of angularly inhomogeneous 
circular cylinder, (a) TM. (b) TE. 
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